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INTRODUCTION 


/  The  use  and  analysis  of  the  signatures  of  light  reflected  by  scattering  sites  in  a 
moving  medium  has  been  actively  pursued  for  the  determination  of  flow  velocities.*'**' 
A  good  overview  of  the  subject  was  given  by  Adrian.*7*'  Here  we  address  the 
problem  of  the  associated  data  reduction,  in  particular  when  the  number  of  scatterers 
is  relatively  small  (the  order  of  100),  and  the  particles  are  relatively  large,  of  the  order 
of  tens  of  microns;  this  is  called '"particle  image  velocimetry  (PIV).^  r  . 

Particle  image  velocimetry  is  of  considerable  interest,  especially  for  flows  where 
there  is  only  limited  transverse  motion,  i.e..  where  particles  in  the  flow  field  generally 
stay  within  well  defined  hyperplanes.  In  many  situations  of  interest  to  ballisticians 
these  conditions  are  met.  Spatial  particle  distributions  within  flows  are  of  some 
importance  since  they  influence  many  processes  of  interest,  including  heat  transfer  and 
velocity  distributions  of  the  carrier  gas  near  bounding  surfaces.  In  addition,  since 
excellent  time  resolution  in  the  data  acquisition  is  now  possible  with  the  use  of  pulsed 
lasers,  the  simultaneous  determination  of  the  position  and  subsequent  inference  of 
velocities  of  a  large  number  of  particles  at  the  same  time  is  a  considerable 
advancement  over  Laser  Doppler  Anemometry  (LDA),5  where  a  large  number  of 
measurements  at  the  same  location  need  to  be  performed  to  establish  the  required  Mow 
statistics  of  a  single  point.  Also,  average  measurements  may  not  accurately  reflect  the 
true  nature  of  the  flow  behavior. 

Going  beyond  current  practice  of  recording  the  data  on  photographic  film 
followed  by  manual  data  reduction,  we  address  the  case  when  two  frames  of  a  scene, 
separated  by  a  small  temporal  displacement,  are  available  from  charge  coupled  device 
(CCD)  recordings.  It  is  shown  that,  even  in  the  case  of  hundreds  of  particles, 
identification,  tracking,  and  determination  of  particle  velocities  can  be  obtained  within 
minutes  even  on  a  slower  mainframe  computer.  In  the  next  section,  a  new  algorithm 
for  PIV  data  reduction  is  described.  This  is  followed  by  an  example  of  a  simulated 
Poiseuille  flow  with  one  hundred  embedded  particles  of  identical  shape.  A  random 
motion  has  been  superimposed  on  each  of  the  particles  as  they  traverse  the  flow  field. 
We  conclude  with  an  assessment  of  the  developed  technique. 


GIVEN  TRACKING  PROBLEM 

The  particles  under  consideration  can  change  in  apparent  shape  and  brightness 
during  their  motion,  since  they  move  in  and  out  of  focus  and  also  might  be  abrading. 
Therefore,  any  potential  tracking  algorithm  cannot  cue  on  individual  particle  signatures 
such  as  edge  segments,  shape,  size  or  optical  brightness.  The  particles  must  be 
considered  identical.  The  main  guides  to  matching  corresponding  particles  must 
therefore  be  their  positions  and  local  configurations.  The  problem  is  further 
exacerbated  by  the  fact  that  only  two  images  are  to  be  given.  Therefore,  a  detailed 
history  of  particle  shape  and/or  trajectory  changes  cannot  be  kept.  We  summarize  the 
problem  as  follows: 

1.  The  particles  undergo  nonrigid  motion  with  a  strong  random  component. 

2.  Many  particles  are  present,  the  order  of  100  or  more.  •  Hence,  the  field  is 
relatively  densely  populated. 


3.  The  particles  are  identical,  and  placed  in  a  uniform  field.  We  equivalently 
model  the  particles  as  uniform  circular  disks  against  a  (different)  uniform 
background. 

4.  Only  two  images  are  given. 

5.  The  time  interval  between  images  is  large  enough  that  particle  travel  of  many 
diameters  can  occur. 

6.  The  deterministic  component  of  motion  is  Poiseuille. 

7.  The  number  of  particles  is  conserved,  over  the  two  images. 

8.  The  aim  is  to  establish  the  particle  correspondences  from  image  1  to  image  2. 
Then  the  difference  in  corresponding  particle  positions  gives  the  required 
velocities. 

One  obvious  method  of  establishing  the  correspondences,  matched  filtering,6  does 
not  apply  because  of  premise  3  above.  The  method  of  optical  flow,6-8  must  also  be 
considered.  It  applies  when  the  image  scene  is  continuously  changing  both  temporally 
and  spatially.  For  example,  a  camera  is  slowly  panning  across  an  image  field. 
However  this  requirement  violates  premises  3  and  5  above,  according  to  which  the 
image  field  is  binary,  and  discontinuous  spatially,  and  the  particle  motion  is  large, 
hence  discontinuous  temporally. 

The  open  literature  does  not  extensively  treat  this  kind  of  tracking  problem. 
Published  work  on  dynamic  scene  analysis. 9-12  generally  limit  attention  to  either  rigid  or 
quasi-rigid  motion  of  a  single,  or  at  most  a  few,  objects.  This  violates  premises  1  and 
2  above.  Other  authors,  e.g.,  as  in  Ref.  13.  allow  for  multiple  particles  but  consider 
cases  where  the  field  is  rather  sparsely  populated,  violating  premise  2.  and  where 
small  particle  travel  occurs  between  frames,  violating  premise  5.  Under  these 
conditions,  there  is  no  correspondence  problem.  Estimation  of  the  particle  positions, 
and  outputing  a  dot  at  each  such  position  on  an  image  screen,  directly  gives  the 
particle  orbits.13  An  ingenious  method  of  tracking  multiple  particles.14  is  based  on  the 
idea  of  predicting  each  particle's  motion  based  on  its  past  history.  However,  this 
violates  premise  4.  according  to  which  there  is  no  extensive  past  history. 

An  analog  method  of  multiple  particle  tracking  was  proposed,15  which  uses  a  far- 
field,  double  exposed  hologram  as  the  recording  medium  for  the  two  images.  However, 
this  approach  is  limited  in  scope  to  two  special  scenarios:  When  the  particles  are 
close  to  uniformly  distributed  over  the  images,  or  the  particle  movements  do  not 
correlate  with  position.  These  are  overly  restrictive  conditions  for  the  given  problem. 

Much  previous  work  in  tracking  has  been  motivated  by  the  needs  of  the  home 
television  industry,  or  by  the  problem  of  inferring  3-D  shape  from  2-D  projections. 
These  have  very  little  in  common  with  the  premises  1  -  8  at  hand.  By  contrast, 
medical  blood  cell  tracking  has  much  in  common  with  our  problem.  The  first 
publication  of  this  kind  appears  to  be  Ref.  16,  where  a  system  for  tracking  single 
white  blood  cells  is  described.  Then,  in  Ref.  17.  an  advance  on  this  work  was 
reported,  whereby  the  automatic  tracking  of  many  ceils  was  carried  out.  The  aim. 
however,  was  not  to  establish  the  particle  correspondences,  but  rather  to  estimate  the 
mean  or  aggregate  motion  of  the  blood  cells.  Hence,  when  multiple  candidates  for  a 
correspondence  occur,  one  is  arbitrarily  selected:  these  events  occur  so  infrequently 
that  the  overall  estimate  of  aggregate  motion  is  not  seriously  affected.  By  contrast, 
because  of  premise  2  above,  such  events  will  frequently  happen  in  our  problem. 
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Also,  each  wrong  decision  on  a  correspondence  gives  directly  a  wrong  velocity 
estimate.  Hence,  we  pay  careful  attention  to  the  problem  of  multiple  candidates. 

Another  medical  tracking  problem  analogous  to  ours  arises  in  2-D  gel 
electrophoresis  imagery.  There,  the  aim  is  to  match  corresponding  protein  spots  from 
one  gel  image  to  another.  The  state-of-the-art  matching  procedure  appears  to  be  that 
of  Ip  and  Potter;1*  see  also  Refs.  19  and  20.  Their  algorithm  consists  of  a  global 
linear  transformation  applied  to  one  image,  followed  by  final,  local  matching  with  the 
second  by  use  of  a  "chamfering"  technique.  The  linear  transformation  takes  advantage 
of  the  powerful  prior  knowledge  that  to  good  extent  the  two  images  are  simply 
translated  and  rotated  versions  of  one  another.  Unfortunately,  this  is  not  true  of 
particles  undergoing  Poiseuille  flow.  Also,  this  linear  step  requires,  for  its  execution, 
the  manual  matching  of  many  reference  spots,  e.g.,  15  to  20  pairs  in  experiments  cited. 
This  conflicts  with  our  aim  of  a  completely  automated  tracking  procedure. 

Skolnick21  has  proposed  an  automatic  matching  approach  in  gel  electrophoresis 
imagery.  Correspondences  are  established  by  a  matching  of  graphs  generated  by 
connecting  each  gel  spot  center  (a  "node")  with  lines  ("edges")  to  its  nearest  neighbors. 
If  the  graphs  for  corresponding  spots  in  successive  frames  are  similar  enough,  and  if 
they  are  not  too  separated  spatially,  then  the  spots  are  accepted  as  matched.  This  was 
a  proposed  technique,  without  demonstration.  Evidently,  it  depends  for  its  utility  on  a 
situation  where  most  particles  preserve  their  relative  positions  during  motion.  This 
violates  our  premises  1  and  6,  according  to  which  the  particle  motion  is  nonuniform. 

Finally.  Greaves22  uses  a  simple  nearest-neighbor  rule,  of  Euclidean  measure,  to 
select  corresponding  microorganisms  in  two  microscopic  slides  viewed  by  a  video 
system.  Unfortunately,  this  simple  selection  rule  will  not  work  in  our  problem,  as  will 
become  evident  from  demonstrations  given  later.  The  field  is  so  crowded  that  the 
correspondences  so  established  will  be  strongly  dependent  upon  the  order  with  which 
the  pairings  are  made,  and  too  often  the  nearest  neighbor  to  a  particle  is  the  wrong 
pairing  for  that  particle. 


PHYSICAL  CONSIDERATIONS 

A  correspondence  algorithm  must  be  sought  which  does  not  rely  on  individual 
object  shape  and/or  brightness  cues  (see  premise  3).  The  only  other  cueing 
information  left  is  possible  knowledge  of  particle  dynamics,  i.e.,  flow  characteristics. 
In  fact,  we  have  such  knowledge.  in  our  problem,  the  particles  overall  follow 
Poiseuilie-like  flow,  which  is  laminar  and  deterministic,  with  a  superimposed  random 
component  of  motion  exemplified  by  local  eddies.  Hence,  in  the  net  there  should  be 
strong  correlation  of  motion,  but  only  over  finite  correlation  distances.  Effectively, 
many  neighboring  particles  move  together,  or  "clump."  Hence,  preference  should  be 
given  to  identifications  that  define  commonly  moving  particles,  or  clumps. 
Furthermore,  the  deterministic  component  of  Poiseuille  flow  is  assumed  to  take 
precedent:  given  no  further  information,  if  a  particle  is  observed  to  belong  to  more 
than  one  clump,  we  presume  that  it  should  be  identified  with  the  larger  clump, 
provided  this  does  not  imply  too  large  a  correlation  distance. 

A  second  piece  of  prior  information  at  hand  is  knowledge  of  a  maximum  possible 
motion  displacement  rmax  over  the  field.  This  can  be  deduced  from  knowledge  of  the 
time  interval  r  between  exposures  and  an  estimate  of  the  maximum  possible  velocity  of 
a  particle.  With  this  information,  potential  correspondences  that  would  require  motions 
greater  than  rmax  can  be  ruled  out.  Furthermore,  if  rmax  <,  N/2.  where  N  is  the 
image  width  (and  height),  a  computational  time  and  core  storage  advantage  results  (see 
below). 
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IMPLEMENTATION  OF  PHYSICAL  CONSTRAINTS:  THE  TRACKING  ALGORITHM 


A  mathematical  operation  that  permits  the  incorporation  of  clumping  and 
maximum  distance  considerations  is  that  of  cross  correlation.  See  the  flow  diagram  in 
Figure  l,  step  (c).  If  image  1  is  cross  correlated  with  image  2,  the  global  maximum  in 
the  output  will  occur  at  a  lag.  or  spatial  vector  displacement  between  the  images, 
where  the  largest  number  of  particles  overlap.  This  defines  the  largest  clump  The 
second-largest  maximum  will  occur  at  a  lag  for  which  the  second- largest  number  of 
particles  overlap;  etc,  for  all  significant  maxima  within  the  largest  possible 
displacement  radius  rmax.  Let  there  be  M  in  all.  The  maxima  were  found  by  simply 
centering  a  3x3  cross  shaped  window  upon  each  output  point.  A  maximum  is  defined 
to  occur  when  the  center  point  is  higher  than  its  four  neighbors  within  the  cross.  The 
M  lag  maxima  describe  the  vector  displacements  of  M  clumps  of  particles  from  image 
I  to  image  2.  In  a  general  scenario,  there  will  be  many  clumps  of  commonly  moving 
particles,  each  moving  by  a  generally  different  amount  and  in  a  different  direction. 
Although  the  cross-correlation  outputs  do  not  identify  which  particles  moved  together  to 
form  each  maximum,  they  at  least  identify  the  lags  in  question.  With  these  known,  a 
different  computer  operation  identifies  the  particles  within  each  clump.  This  is 
described  in  the  third  paragraph  following. 


The  cross  correlation  operation  was  also  chosen  because  it  is  convenient  to 
implement  digitally.  Let  image  I  be  denoted  as  f.  image  2  denoted  as  g.  The  Fast 
Fourier  Transform  (FFT)  algorithm  permits  the  cross  correlation  f  e  g  to  be  computed 
for  two  images  f  and  g.  as23 

f  »  g  »  FFT-1  {FFT(0  FFT(g)*}  .  (1) 

The  asterisk  denotes  a  complex  conjugate,  and  the  -1  indicates  an  inverse  FFT.  The 
FFT  algorithm  is  notably  fast.  Also,  by  Eq.  (1).  only  three  FFT  operations  need  be 
done  to  produce  the  required  output  that  locates  the  lag  maxima.  Finally,  size 
considerations  should  be  mentioned.  If  the  images  f  and  g  are  NxN  pixels,  then 

ordinarily  the  output  fug  will  be  2Nx2N  pixels,  i.e..  four  times  the  a^ea  of  each 

input.  However,  if  it  is  known  that  no  particle  has  been  translated  from  f  to  g  by 

greater  than  N/2  pixels,  i.e..  rmax  ^  N/2.  then  all  output  maxima  must  instead  lie 

within  an  NxN  central  field  of  the  output.  This  permits  a  reduction  in  required 
computer  core  storage  by  a  factor  of  8  (factors  of  2  in  each  direction,  plus  a  factor  of 
2  for  the  complex  arithmetic  used). 

A  further  size  reduction  was  enforced.  By  extensive  use  of  disk  intermediary 

storage  files,  it  is  possible  to  require  but  one  core  storage  array  of  dimension 

(2,  N,  N).  Without  the  use  of  such  disk  files,  three  such  arrays  would  have  been 

needed,  plus  seven  of  dimension  (N,  N).  For  an  image  of  size  N  =  128  or  larger,  this 

amounts  to  a  considerable  reduction  in  core  storage  requirement. 

Assume  that,  in  a  pre-processing  step,  every  particle  in  the  two  images  is 
replaced  by  an  identical  disk  of  known  (x,  y)  position.  See  steps  (a),  (b)  in  Figure  1. 
The  particles  within  each  clump  were  then  identified  in  the  following  way.  (See  step 
(e)  of  Figure  1.)  Displace  image  2  from  image  I  by  one  of  the  M  lags  defined  abow. 
Then  for  each  particle  position  in  image  1,  every  particle  position  in  image  2  is 
sampled  to  see  if  the  two  particle  disks  overlap.  Of  all  such  overlapped  disks,  the 
pair  with  the  smallest  mutual  displacement  is  chosen  as  defining  a  particle-pair 
correspondence.  With  the  given  scenario  of  little  a  priori  knowledge  governing  particle 
motion,  it  is  logical  to  prefer  nearest-neighbor  pair  identifications  over  all  candidates 
that  are  otherwise  equally  valid. 


M. 
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(a) 


(b) 


(c) 


(d) 


(e) 


(0 


(g) 


Figure  1.  Flow  diagram  of  tracking  algorithm. 


The  M  lags  are  ordered  such  that  their  corresponding  cross  correlation  maxima 
increase  in  size,  before  the  foregoing  identifications  are  made.  See  step  (d)  in  Figu:e 
1.  Then  at  each  lag,  a  new  set  of  particle  correspondences  is  established  as  abo' e. 
with  each  such  set  replacing  any  conflicting  correspondences  established  at  previous 
lags.  For  example,  if  at  lag  1  particle  10  of  image  1  is  identified  with  particle  8  of 
image  2.  i.e..  identification  (10,  8)  is  made;  while  at  lag  2  identification  (10.  5)  is 
made;  then  identification  (10.  5)  takes  precedence,  until  possibly  at  some  subsequent  lag 
particle  10  is  assigned  to  yet  a  different  particle  of  image  2.  In  this  way,  the  larger 
clumps  associated  with  the  later  lags  are  given  preference  in  defining  particle 
identifications. 
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As  with  the  nearest-neighbor  rule  described  above,  this  is  intuitively  a  correct 
choice.  If  the  particle  flow  was  perfectly  smooth  and  laminar,  all  particles  would 
travel  together  in  one  grand  clump.  However,  the  extent  to  which  random  eddies  and 
other  sources  of  randomness  actually  creep  in  is  not  known.  Hence,  it  makes  sense  to 
give  preference  to  large  scale,  or  laminar,  flow. 

Notice  that  in  the  previous  step  particle  8  in  image  2  was  "bumped"  as  a 

possible  candidate  for  an  identification.  After  looping  over  all  lags  in  the  previous 
step,  there  will  be  a  list  of  such  bumped  particles  from  image  2.  Likewise,  there  will 
be  a  list  of  still-unassigned  particles  from  image  1  (those  that  did  not  have  an 
overlapped  disk,  for  any  lag,  with  a  particle  in  image  2).  These  two  lists  must 

somehow  be  matched  up  1:1.  At  this  point,  all  clumping  tendencies  have  been 

established.  The  only  rule  left  for  match  up  is  the  nearest-neighbor  rule.  Hence,  we 
loop  over  all  such  image  1  particles,  finding  the  closest  image  2  particle,  and  then 
eliminating  each  from  its  list.  (See  step  (~)  of  Figure  1).  This  is  admittedly  an 
imperfect  procedure,  since  the  identifications  will  depend  somewhat  upon  the  order  in 
which  the  particles  are  processed.  For  example,  the  particle  pair  (3,  8)  may  define  the 
smallest  particle  distance  to  particle  3  of  image  1,  but  8  may  have  already  been 
assigned  to  particle  2  of  image  1  on  the  basis  of  closeness  to  2.  In  practice,  this  kind 
of  error  does  not  happen  very  often.  However,  to  minimize  its  occurrence  a 
subsequent  processing  step  is  taken,  (g)  of  Figure  1. 

In  this  step,  each  identification  (m|.  n().  i  -  1 . P.  where  P  is  the  total 

number  of  particles  in  each  image,  is  compared  with  every  other  identification  (m  ,  n^), 
j  #  i.  to  see  if  an  interchange  of  identifications,  to  (mj.  nj).  (mj,  n()  will  produce  a 
smaller  total  displacement  distance  over  both  pairings.  If  it  does,  then  the  interchange 
is  made.  Again,  this  rule  gives  preference  to  small  displacements.  In  practice,  the 
rule  eliminates  all  or  nearly  all  errors  described  in  the  previous  paragraph. 


DEMONSTRATION 

Figure  2  shows  a  simulated  case  study  of  particle  identification.  To  aid 
visualization,  image  1  consists  of  a  regular  10x10  grid  of  particles,  shown  as  white 
diamonds  (discrete  versions  of  circles).  Field  size  is  128x128  pixels.  The  100  particles 
are  simulated  to  obey  randomized  Poiseuille  flow,  generally  to  the  right.  The  center 
row  of  particles  move  maximally  and  the  top  and  bottom  rows  minimally,  according  to 
a  parabolic  dependence  on  row.  (The  top  and  bottom  rows  are  taken  to  be  near  the 
walls  of  a  pipe  containing  the  flow.  Because  of  frictional  interactions  with  the  walls, 
particles  close  to  them  tend  to  move  less  than  those  in  the  main  flow  in  the  middle.) 
The  center  row  has  the  maximum  parabolic  displacement,  value  rmax.  Parameter  rmax 
was  given  the  value  6.0  pixels. 

Each  particle's  direction  of  flow  is  made  to  depart  randomly  from  the  horizontal, 
with  maximum  randomness  again  at  (he  walls  since  this  is  where  eddy  currents  tend 
to  maximally  occur.  The  result  is  the  black  particle  positions  in  Figure  2.  As 
examples  of  extreme  directional  randomness,  notice  the  lower-left  most  black  particle 
positions.  The  black  particles  comprise  image  2. 

These  two  images  were  then  processed,  according  to  tracking  algorithm  (a) 
through  (g)  of  Figure  1.  At  step  (c),  M  =  5  lag  maxima  were  found  within  the 
feasible  region  defined  by  rmax.  The  result  is  the  correspondences  shown  by  black 
connecting  lines  in  Figure  2.  (Where  a  black  pixel  and  a  while  one  overlap,  as  in  the 
upper-left  most  pair,  the  black  one  dominates.)  These  results  are  encouraging,  since 
every  particle  was  correctly  tracked. 
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Figure  2.  Simulated  tracking  problem.  White  particles  are  in  image  1, 
black  particles  in  image  2.  Correspondences  established  by  the  tracking 
algorithm  are  designated  by  black  connecting  lines. 

It  is  interesting  to  compare  these  results  with  what  a  simple  nearest-neighbor  rule 
might  have  given.  On  the  basis  of  nearest  neighbors,  some  erroneous  identifications 
would  be  made,  and  some  would  be  ambiguous.  For  example,  consider  the  white 
particle  in  row  6,  column  9.  Its  nearest  black  neighbor  is  to  its  left.  6  pixels  away. 
However,  the  tracking  algorithm  correctly  paired  it  ir^tead  with  the  particle  to  its 
right,  distance  (62  +  32)1/2  >  6.3  pixels  away.  Although  this  distance  is  larger,  it  was 
identified  by  the  program  with  a  large  clump  of  commonly  moving  particles,  and  hence 
was  given  preference  at  step  (e)  of  the  algorithm.  Also,  for  many  particles  nearest- 
neighbors  would  have  led  to  an  impasse,  since  two  white  particles  exist  which  are  the 
same  distance  from  a  given  black  one.  For  example,  consider  the  white  particle 
located  at  row  3,  column  3.  The  black  particles  to  its  immediate  left  and  right  are 
each  distance  (62  +  I2)1/2  from  it.  A  supplemental y  rule  would  have  had  to  be 
invented  in  order  to  resolve  the  ambiguity. 


DISCUSSION 

Although  the  results  in  Figure  2  are  typical,  the  algorithm  does  not  of  course 
always  work  this  well.  The  central  consideration  is  the  amount  of  particle  motion 
compared  with  the  interparticle  distances.  If  a  black  particle  has  moved  so  much  that 
it  overtakes  the  white  particle  immediately  to  its  right,  it  might  be  erroneously 
identified  with  the  latter  white  particle.  Hence,  particle  movement  should  be  less  than 
the  interparticle  distances  in  image  I. 


Another  factor  of  importance  is  the  designated  size  of  the  disks.  Notice  that  disk 
size,  as  employed  in  the  algorithm,  does  not  have  to  be  the  actual  or  physical  disk 
size.  In  practice,  disk  size  only  enters  into  the  algorithm  in  step  (e),  where  it  is  used 
as  a  kind  of  interaction  distance.  Specifically,  it  defines  candidate  correspondences  for 
an  image  1  particle.  To  be  a  candidate,  an  image  2  particle  must  overlap  the  disk  of 
the  image  1  particle.  (The  closest  such  image  2  particle  is  then  selected.)  Obviously, 
then,  the  size  taken  for  the  disk  governs  the  chosen  set  of  candidates.  With  too  small 
a  disk  size,  such  a  small  set  of  candidates  might  be  generated  that  the  correct 
candidate  is  missed.  Or,  too  large  a  disk  will  too  strongly  favor  the  last  lags  used  in 
step  (e),  i.e.  will  tend  to  make  correspondences  that  all  correspond  to  maximum  clump 
size.  Hence,  disk  size  is  a  useful  tuning  parameter.  A  good  disk  size  is  identified  by 
reasonable  looking  correspondences,  in  the  judgment  of  the  user.  For  the  results 
shown  in  Figure  2,  the  chosen  disk  radius  was  2  pixels,  as  shown.  This  was  the 
second  run  of  the  problem.  In  the  first,  a  disk  radius  of  1  pixel  was  used  instead. 
This  resulted  in  all  correspondences  but  eight  in  row  9  being  correct,  a  success  rate  of 
92/100.  Hence,  the  algorithm  is  fairly  tolerant  of  disk  size.  A  range  of  disk  sizes 
will  give  good  results. 

The  simulation  shown  in  Figure  2  was  carried  out  by  a  Cyber  175  mainframe 
computer.  CPU  time  was  51  sec,  and  central  memory  required  60,200  octal  words  in 
total  (including  all  program  statements  and  all  execution  arrays.)  A  major  part  of  the 
CPU  time  was  taken  by  the  seven  FFT  operations  performed— three  to  locate  particle 
positions  in  image  1,  three  similarly  for  image  2.  and  one  more  to  cross  correlate 
image  1  with  image  2.  For  128x128  pixel  images,  it  was  empirically  found  that  CPU 
time  t  varies  with  particie  number  P  according  to  a  linear  relation 

t  -  P/3  +  20  sec.  1  £  P  ^  100.  (2) 

Since  the  FFT  operations  are,  of  course,  independent  of  P.  the  20  sec  contribution  to 
Eq.  (2)  must  be  due  to  them.  All  other  operations  in  the  tracking  algorithm  require 
P/3  sec.  over  the  indicated  range  of  P.  This  bodes  well  for  applications  to  larger 
problems,  where  P  might  be  the  order  of  1000  particles,  providing  that  approximate 
linearity  holds  for  these  P  values  as  well. 


CONCLUSIONS 

A  tracking  algorithm  has  been  developed  that  satisfactorily  tracks  100  or  more 
identical  particles.  Time  and  storage  requirements  are  modest,  with  CPU  time 
approximately  linear  in  the  number  of  particles.  For  proper  use.  particle  movements 
should  not  exceed  interparticle  distances  in  image  1.  Program  operation  is  adjustable 
by  an  input  effective  particle  size  parameter,  which  is  varied  until  satisfactory 
correspondences  are  made.  Application  of  the  algorithm  to  1000  or  more  particles 
seems  a  reasonable  future  prospect. 
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ATTN:  ATZK-CD-M5 
M.  Falkovitch 
Armor  Agency 
Fort  Knox,  KY  40121-5215 

2  Commander 

Naval  Sea  Systems  Command 
ATTN:  SEA62R 
SEA  64 

Washington,  DC  20362-5101 

1  Commander 

Naval  Air  Systems  Command 
ATTN:  AlR-954-Tech  Lib 
Washington,  DC  20360 

1  Assistant  Secretary  of  the 

Navy  (R,  E,  and  S) 

ATTN:  R.  Reichenbach 
Room  5E787 
Pentagon  Bldg. 

Washington,  DC  20350 

1  Naval  Research  Lab 
Tech  Library 
Washington,  DC  20375 

5  Commander 

Naval  Research  Laboratory 
ATTN:  L.  Harvey 

J.  McDonald 
E.  Oran 
J.  Shnur 
R.  Doyle 

Washington,  D.C.  20375 

5  Commander 

Naval  Surface  Weapons  Center 
ATTN:  Code  G33,  J.  L.  East 
W.  Burrell 
J.  Johndrow 
Code  G23,  D.  McClure 
Code  DX-21  Tech  Lib 
Dahlgren,  VA  22448-5000 

2  Comander 

US  Naval  Surface  Weapons 
Center 

ATTN:  J.  P.  Consaga 
C.  Gotzmer 

Indian  Head,  MD  20640-5000 


4  Commander 

Naval  Surface  Weapons  Center 
ATTN:  S.  Jacobs/Code  240 
Code  730 

K.  Kim/Code  R-13 
R.  Bernecker 

Silver  Spring,  MD  20903-5000 

2  Commanding  Officer 

Naval  Underwater  Systems 
Center 

Energy  Conversion  Dept. 

ATTN:  CODE  5B331,  R.  S.  Lazar 
Tech  Lib 

Newport,  RI  02840 

4  Commander 

Naval  Weapons  Center 
ATTN:  Code  388,  R.  L.  Derr 
C.  F.  Price 
T.  Boggs 
Info.  Sci.  Div. 

China  Lake,  CA  93555-6001 

2  Superintendent 

Naval  Postgraduate  School 
Dept,  of  Mech.  Engineering 
Monterey,  CA  93943-5100 

1  Program  Manager 

AFOSR 

Directorate  of  Aerospace 
Sciences 

ATTN:  L.  H.  Caveny 
Bolling  AFB,  DC  20332-0001 

5  Commander 

Naval  Ordnance  Station 
ATTN:  P.  L.  Stang 
L  Torrcyson 
T.  C.  Smith 
D.  Brooks 
Tech  Library 

Indian  Head,  MD  20640-5000 

I  AFSC/SDOA 

Andrews  AFB,  MD  20334 

3  AFRPL/DY,  Stop  24 

ATTN:  J.  Levine/DYCR 

R.  Corlcy/DYC 
D.  Williams/DYCC 
Edwards  AFB,  CA  93523-5000 

1  AFRPL/TSTL  (Tech  Library) 

Stop  24 

Edwards  AFB,  CA  93523-5000 

I  AFATL/DLYV 

Egiin  AFB,  FL  32542-5000 
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1 


I 


I 


1 


AFATL/DLXP 

Eglin  AFB,  FL  32542-5000 


2 


AFATL/DLJE 

Eglin  AFB,  FL  32542-5000 


AFATL/DOIL 

ATTN:  (Tech  Info  Center) 

Eglin  AFB,  FL  32542-5438 


1 


NASA 

Langley  Research  Center 
ATTN:  G.B.  Northam/MS  168 
Langley  Station 
Hampton,  CA  23365 

NASA/Lyndon  B.  Johnson  Space 
Center 

ATTN:  NHS-22,  Library 
Section 

Houston,  TX  77054 


1 


1 


Calspan  Corporation 
ATTN:  C.  Morphy 
P.  O.  Box  400 
BufTalo,  NY  14225-0400 

General  Electric  Company 
Armament  Systems  Dept. 
ATTN:  M.  J.  Bulman, 

Room  131 1 

128  Lakeside  Avenue 

Burlington,  VT  05401-4985 

General  Applied  Sciences  Lab 
ATTN:  J.  Erdos 
Merrick  &  Stewart  Avenues 
Westbury  Long  Isld,  NY  1 1 590 

General  Motors  Rsch  Labs 
Physics  Department 
ATTN:  R.  Teets 
Warren,  Ml  48090 


1  AFELM,  The  Rand  Corporation 
ATTN:  Library  D  (Required  or 
1700  Main  Street  Classified 
Santa  Monica  CA  Only) 
90401-3297 

2  AAI  Corporation 
ATTN:  J.  Hebert 

J.  Frankie 
P.  O.  Box  6767 
Baltimore,  MD  21204 


1  11TRI 

ATTN:  M  J.  Klein 
10  W.  35th  Street 
Chicago,  IL  60616-3799 

1  Hercules  Inc. 

Allegheny  Ballistics 
Laboratory 
ATTN:  R.  B.  Miller 
P.  O.  Box  210 

Cumberland,  MD  21501  -02 1 0 


1 


1 


1 


1 


1 


Aerojet  Ordnance  Company 
ATTN:  D.  Thatcher 
2521  Michelle  Drive 
Tustin,  CA  92680-7014 

Aerojet  Solid  Propulsion  Co. 
ATTN:  P.  Michcli 
Sacramento,  C A  95813 


Atlantic  Research  Corporation 
ATTN:  K.  King 
5390  Cheorokee  Avenue 
Alexandria,  V A  22312-2302 


AVCO  Everett  Rsch  Lab 
ATTN:  D.  Stickler 
2385  Revere  Beach  Parkway 
Everett,  MA  02149-5936 


Battelle  Memorial  Institute 
ATTN:  Tech  Library 
505  King  Avenue 
Columbus,  OH  43201-2693 


1  Hercules,  Inc. 

Bacchus  Works 
ATTN:  K.  P.  McCarty 
P.  O.  Box  98 
Magna,  UT  84044-0098 

1  Hercules,  Inc. 

Radford  Army  Ammunition  Plant 
ATTN:  J.  Pierce 
Radford,  V  A  24141-0299 

1  Honeywell,  Inc.  -  MN64  2200 

Defense  Systems  Division 
ATTN:  C.  Hargreaves 
6110  Blue  Circle  Drive 
Minnetonka  MN  55436 

1  Institute  of  Gas  Technology 

ATTN:  D.  Gidaspow 
3424  S.  State  Street 
Chicago,  IL  60616-3896 

l  Lawrence  Livermore  National 

Laboratory 
ATTN:  C.  Westbrook 
P.O.  Box  808 
Livermore,  CA  94550 
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1  Lawrence  Livermore  National 

Laboratory 
ATTN:  L-355, 

A.  Buckingham 
M  Finger 
P.  O.  Box  808 

Livermore,  CA  94550-0622 

1  Lawrence  Livermore  National 

Laboratory 

ATTN:  L-324/M  Constantino 
P.  O.  Box  808 

Livermore,  CA  94550-0622 

3  Los  Alamos  National  Lab 

ATTN:  B.  Nichols 

T7,  MS-B284 
C.  Mader 
K.  Hanson 
P.O.  Box  1663 
Los  Alamos,  NM  87545 

1  Massachusetts  Institute  of 
Technology 

Dept  of  Mechanical 
Engineering 
ATTN:  T.  Toong 
77  Massachetts  Avenue 
Cambridge,  MA  02139-4307 

2  Mathematics  Research  Center 
ATTN:  B.  Noble 

J.  Nohcl 
610  Walnut  Street 
Madison,  WI  53706 

1  Olin  Corporation 

Badger  Army  Ammunition  Plant 
ATTN:  R.  J.  Thiede 
Baraboo,  Wl  53913 

1  Olin  Corporation 

Smokeless  Powder  Operations 
ATTN:  D.  C.  Mann 
P.O.  Box  222 

St.  Marks,  FL  32355-0222 

1  Paul  Gough  Associates,  Inc. 

ATTN:  P.  S.  Gough 
P.  O.  Box  1614, 

1048  South  St. 

Portsmouth,  NH  03801-1614 

1  Princeton  Combustion  Research 

Lab.,  Inc. 

ATTN:  N.A.  Messina 
475  US  Highway  One 
Monmouth  Junction,  NJ 
08852-9650 


No.  of 

Copies  Organization 

1  Science  Applications,  Inc. 

ATTN:  R.  B.  Edelman 
23146  Cumorah  Crest  Drive 
Woodland  Hills,  CA  91364-3710 

3  Thiokol  Corporation 

Huntsville  Division 
ATTN:  D.  Flanigan 
R.  Glick 
Tech  Library 
Huntsville,  AL  35807 

2  Thiokol  Corporation 
Elkton  Division 
ATTN:  R.  Biddle 
Tech  Lib. 

P.  O.  Box  241 
Elkton,  MD  21921-0241 

2  United  Technologies 

Chemical  Systems  Division 
ATTN:  R.  Brown 
Tech  Library 
P.  O.  Box  358 

Sunnyvale,  CA  94086-9998 

1  Veritay  Technology,  Inc. 

ATTN:  E.  Fisher 
4845  Millersport  Hwy. 

P.  O.  Box  305 

East  Amherst,  NY  14051-0305 

1  Universal  Propulsion  Company 

ATTN.  H.  J.  McSpaddcn 
Black  Canyon  Stage  1 
Box  1140 

Phoenix,  AZ  85029 

1  Brigham  Young  University 

Dept,  of  Chemical  Engineering 
ATTN:  M  Beckstead 
Provo,  UT  84601 

1  California  Institute  of  Tech 

204  {Carman  Lab 
Main  Stop  301-46 
ATTN:  F.  E.  C.  Culick 
1201  E.  California  Street 
Pasadena,  C A  91109 

1  California  Institute  of  Tech 

Jet  Propulsion  Laboratory 
ATTN:  L.  D.  Strand 
4800  Oak  Grove  Drive 
Pasadena,  C  A  911 09-8099 

1  University  of  California, 

Berkeley 

Mechanical  Engineering  Dept. 
ATTN:  J.  Daily 
Berkeley,  CA  94720 
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1  University  of  Florida 

Dept,  of  Chemistry 
ATTN:  J.  Winefordner 
Gainesville,  FL  32611 

1  University  of  Illinois 

Dept  of  Mech/Indust  Engr 
ATTN:  H.  Krier 
144  MEB;  1206  N.  Green  St. 
Urbana, IL  61801-2978 

1  University  of  Massachusetts 

Dept,  of  Mech.  Engineering 
ATTN:  K.  Jakus 
Amherst,  MA  01002-0014 

1  University  of  Minnesota 

Dept,  of  Mech.  Engineering 
ATTN:  E.  Fletcher 
Minneapolis,  MN  55414-3368 

1  Case  Western  Reserve 

University 

Division  of  Aerospace 
Sciences 
ATTN:  J.  Tien 
Cleveland,  OH  44135 

3  Georgia  Institute  of  Tech 

School  of  Aerospace  Eng. 
ATTN:  B.  T.  Zmn 
E.  Price 
W.  C.  Strahle 
Atlanta,  GA  30332 

1  Virginia  Polytechnic 

Institute  and 
State  University 
ATTN:  J.A.  Schctz 
R.  T.  Smith 

Blacksburg,  VA  24061 

1  Institute  of  Gas  Technology 

ATTN:  D.  Gidaspow 
3424  S.  State  Street 
Chicago,  IL  60616-3896 

1  G.  M.  Faeth 

Pennsylvania  State  University 
Applied  Research  Laboratory 
University  Park,  PA 
16802-7501 

1  Pennsylvania  State  University 

Dept,  of  Mech.  Engineering 
ATTN:  K.  Kuo 
University  Park,  PA 
16802-7501 


No.  of 

Copies  Organization 

1  Princeton  University 
MAE  Dept 

ATTN:  F.  A.  Williams 
Princeton,  NJ  08544 

2  Purdue  University 
School  of  Mechanical 

Engineering 

ATTN:  N.  M.  Laurendeau 
S.  N.  B.  Murthy 
TSPC  Chaffee  Hall 
West  Lafayette,  IN  47907-1199 

I  Purdue  University 

School  of  Aeronautics 
and  Astronautics 
ATTN:  J.  R.  Osborn 
Grissom  Hall 
West  Lafayette,  IN  47906 

1  Rensselaer  Polytechnic  Inst. 
Department  of  Mathematics 
lroy,  NY  12181 

2  Director 

Los  Alamos  Scientific  Lab 
ATTN:  T3,  D.  Butler 

M  Division,  B.  Craig 
P.O.Box  1663 
Los  Alamos,  NM  87544 

1  Stevens  Institute  of 

Technology 
Davidson  Laboratory 
ATTN:  R.  McAlevy,  III 
Castle  Point  Station 
Hoboken,  NJ  07030-5907 

1  Rutgers  University 

Dept,  of  Mechanical  and 
Aerospace  Engineering 
ATTN:  S.  Temkin 
University  Heights  Campus 
New  Brunswick,  NJ  08903 

1  University  of  Southern 

California 

Mechanical  Engineering  Dept. 
ATTN:  0HE200,  M.  Gerstcin 
Los  Angeles,  CA  90089-5199 

2  University  of  Utah 

Dept,  of  Chemical  Engineering 
ATTN:  A.  Baer 

G.  Flandro 

Salt  Lake  City,  UT  84112-1194 

1  Washington  State  University 

Dept,  of  Mech.  Engineering 
ATTN:  C.  T.  Crowe 
Pullman,  WA  99163-5201 
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Dir,  USAMSAA 

ATTN:  AMXSY-D 

AMXSY-MP,  H.  Cohen 

Cdr,  USATECOM 

ATTN:  AMSTE-SI-F 

AMSTE-CM-F,  L.  Nealley 


Cdr,  CSTA 

ATTN:  STECS-AS-H,  R.  Hendricksen 


Cdr,  CRDC,  AMCCOM 

ATTN:  SMCCR-RSP-A 


SMCCR-MU 

SMCCR-SPS-IL 


Organization 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the 
reports  it  publishes.  Your  comment s/answers  to  the  itews/questions  below  will 
aid  us  in  our  efforts. 

1.  BRL  Report  Number _ Date  of  Report _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or 

other  area  of  interest  for  which  the  report  will  be  used.) _ 


4.  How  specifically,  is  the  report  being  used?  (Information  source,  design 
data,  procedure,  source  of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 
as  man-hours  or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved, 
etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 


Name 


CURRENT 

ADDRESS 


Organization 


Address 


City,  State,  Zip 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


Name 


Organization 

ADDRESS  _ 

Address 


City,  State,  Zip 


(Remove  this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  mail.) 


FOLD  HERE 


Director 

U.S.  Army  Ballistic  Research  Laboratory 
ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-5066 


OFFICIAL  BUSINESS 

PENALTY  PON  PRIVATE  USE.  SSOO 


BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  NO  12062  WASHINGTON, DC 
POSTAGE  WILL  BE  PAID  BY  DEPARTMENT  Of  THE  ARMY 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


Director 

U.S.  Army  Ballistic  Research  Laboratory 
ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-9989 


FOLD  HERE 


